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Abstract. We present a new multi-fluid, grid MHD code PIERNIK, 
which is based on the Relaxing TVD scheme ( |Jin and Xin, 1995| . The 
original scheme (see Trac & Pen (|20U3p and Pen et al. (|2003p ) has been 
extended by an addition of dynamically independent, but interacting 
fluids: dust and a diffusive cosmic ray gas, described within the fluid 
approximation, with an option to add other fluids in an easy way. 
The code has been equipped with shearing-box boundary conditions, 
and a selfgravity module, Ohmic resistivity module, as well as other 
facilities which are useful in astrophysical fluid-dynamical simulations. 
The code is parallelized by means of the MPI library. In this paper we 
present an extension of PIERNIK, which is designed for simulations 
of diffusive propagation of the Cosmic-Ray (CR) component in the 
magnetized ISM. 



1 Cosmic Ray transport 

The CR-MHD extension of PIERNIK code (Hanasz et al. (|2009ap . (|2009bp . 
(|2009cp ) is aimed at studies of the magnetohydrodynamical dynamo process in- 
duced by buoyancy of CRs in stratified atmospheres of galactic disks (Parker 
(|1992p . Hanasz et al. (|2004p ') investigated previously, in the shearing-box ap- 
proximation, with the aid of ZEUS-3D code, extended with the CR transport 
algorithm (Hanasz & Lesch (|2003p ). 

To describe cosmic-ray (CR) propagation in the interstellar medium (ISM) we 
use the diffusion-advection equation (see Schlickeiser & Lerche (|1985p ) 

dte + V{ev) = -pV ■ v + V{kVe) + Qcr (1-1) 

together with the adiabatic equation of state for cosmic rays 

Per = (7cr - l)ecr, (1-2) 
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in addition to the standard set of MHD equations (Hanasz & Lesch (|2003p V The 
source term Qcr on the rhs. of Eqn. p.ip corresponds to the production of CRs 
in supernova remnants. The diffusion term is written in the tensorial form to ac- 
count for anisotropic diffusivity of CRs, where K is the diffusion tensor describing 
magnetic field-aligned CR diffusion (see Ryu et al (|2003p ) 

Kij = K^5ij + (i^ii - Ki_)ntnj, = Bt/B, (1.3) 

We note that in the presence of CRs an additional source term: — VPcr should be 
included in the gas equation of motion (see e.g. Berezinski et al. (|1990p ). in order 
to incorporate the effects of CRs on gas dynamics. 

In order to adopt the CR transport equation to the conservative scheme of 
PIERNIK code, we write Eqn. ()l.ip in the conservative form 

dtecr + V • F„,adv + V • Fcr,diff = -PcrV • V + Q„, (1.4) 

where ecr is CR energy density, .Fcr.adv — £ci-v, is the flux of CRs advected by the 
gas flow, -Fcr,diff = — -ftTVccr is CR diffusion flux and Qa is the CR source term. 

The left hand side of the CR transport equation is treated in a conservative 
manner, while the terms on r.h.s. are added as source terms. The advection and 
source steps (pci-V • v) for CRs are implemented within the RTVD scheme, while 
the CR diffusion step is realized outside the Relaxing TVD routine. The update 
of CR energy, corresponding to the diffusion term is performed with the aid of a 
directionally split, explicit algorithm (Hanasz & Lesch (|2003p ). which is first order 
in time and space. The source step corresponding to the injection of CRs in SN 
remnants is realized once per double timestep, outside the directional sweeps of 
fluid updates. 

The explicit CR diffusion algorithm implemented in PIERNIK code is subject 
to the timestep limitation resulting from the von Neumann stability analysis. The 
timestep for the diffusive part of the CR diffusion-advection equation imposed in 
the code is 

„ min(Ax, Aw, Az)^ 

= 1 ^ (1-5) 

A II + Aj_ 

where Ccr < 1 is the Courant number for the CR diffusive transport algorithm. 
Ax, Ay and Az are cell sizes. The numerical stability of the overall CRMHD algo- 
rithm, is achieved by a proper monotonic interpolation of CR gradient components 
computed on cell boundaries (see Hanasz & Lesch (|2003p ). We note, that the ap- 
propriate choice of boundary conditions for the highly diffusive CR component is 
to set a fixed value (zero) of CR energy density on external domain boundaries. 



2 Test problems for CR transport 

To test the magnetic field-aligned cosmic ray diffusion we present a simple 2D 
setup with uniform and diagonal magnetic field in the doubly-periodic computa- 
tional domain. Parameters of the initial setup for the simulation are (in arbitrary 
units): po = 1, po = 1, B^ = 3, By = 3, 7 = 5/3, j„ = 4/3, Xmin = -500, 
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Fig. 1. Diffusion of cosmic rays along an inclined magnetic field: the initial spheroidal 
distribution of ecr at f = and the ellipsoidal distribution at t = 20 and t = 60. The last 
panel shows thermal gas density and velocity vectors at t = 60. The apparent flow of gas 
along the magnetic fleld direction is due to the CR pressure gradient, pushing gas along 
magnetic field lines. A magnetosonic wave, propagating in the direction perpendicular 
to magnetic field is also present. 



a^max = 500, j/niin = —500 and ?/max = 500. At i = a portion of CRs, forming a 
2D Gaussian profile, with half-width equal to 50 units and = 8 at maximum 
around the domain center. The diffusion coefficients are K\\ = 1000 and A'^ = 0. 

The results of the test run demonstrate that the CR diffusion proceeds along 
magnetic field lines, as expected (see first three panels of Fig. [1]). A detailed 
quantitative analysis ensures that in case of passive CR propagation (without the 
back-reaction of CR pressure on the thermal gas) numerical results fit accurately to 
the analytical solution. In the present case of active CR propagation, CR pressure 
gradients affect thermal gas (see the fourth panel of Fig. [1]). The excess of cosmic 
ray pressure near the center of computational domain accelerates gas, along the 
oblique magnetic field, forming an ellipsoidal cavity in the gas distribution. 

The present implementation of CR transport within the very efficient and flexi- 
ble Relaxing TVD scheme (Pen ct al. p003|) ). combined with the MPI parallcliza- 
tion of PIERNIK; makes it possible to study the dynamic of CRs, and CR-drivcn 
dynamo in global simulations of galactic disks. First results of global galactic 
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dynamo simulations (Hanasz et al. (|2009dp . (120096^ ) demonstrate that magnetic 
fields can be efficiently amplified to equipartition values, on the timescale of galac- 
tic rotation, starting from weak magnetic fields of stellar origin. 

In a more general case the cosmic ray (CR) component can be considered as 
an additional set of fluids extending the vector u of conservative variables. A 
subsequent extension of the CR transport module in PIERNIK code, aiming at 
energy dependent treatment of CR-electrons. and incorporation of synchrotron 
losses is currently under development. 
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